Aquest és un exercici de mineria de dades sobre l’estudi dels terratrèmols o activitat sísmica global, és a dir plantegem l’estudi de l’activitat sísmica mundial.
Disposar d’un bon conjunt de dades i un bon coneixement del camp del que volem fer l’estudi ajuda a fer-se bones preguntes i a orientar-se cap a l’objectiu on es vol arribar. En aquest cas, no disposo d’un bon coneixement de geologia i les dades…ja veurem, però sí que personalment, crec és un tema interessant i val la pena dedicar-li aquesta pràctica.
Comencem doncs amb un joc de dades procedents de Kaggle, en particular:
Aquest dataset que hem triat, Global Significant Earthquake Database from 2150BC, és un llistat d’uns 6000 terratrèmols que van del 2150 BC fins a l’actualitat. Això ja ens diu de bones a primeres que les dades prèvies a l’ús d’instruments i recull de dades, probablement són incompletes. De fet, ja avanço que aquest conjunt de dades és altament incomplet.
Com a objectiu general, volem mostrar quines són les zones i països amb més activitat sísmica i perill de tsunamis i el seu impacte en aquestes àrees, això posarà de manifest l’estructura, posició i dinàmica de les plaques tectòniques i a la vegada les possibles conseqüències per la població. Volem veure quins països són els més actius sísmicament. També volem fer una classificació i veure què fa que un terratrèmol sigui devastador, (significant). Intentarem esbrinar relacions o correlacions entre magnitud del sisme, profunditat a la que es produeix i altres variables. El valor d’aquestes dades pot ajudar a minimitzar conseqüències o a preveure i millorar infraestructures. De totes maneres com s’ha vist sempre, les àrees amb menys recursos i més pobres són les més damnificades.
library(tidyverse)
El primer que fem és carregar les dades i fer una ullada per veure de quantes variable, de quin tipus són i de quants registres disposem.
dades <- read_csv('Worldwide-Earthquake-database.csv')
## Rows: 6193 Columns: 47
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): FLAG_TSUNAMI, COUNTRY, STATE, LOCATION_NAME
## dbl (43): I_D, YEAR, MONTH, DAY, HOUR, MINUTE, SECOND, FOCAL_DEPTH, EQ_PRIMA...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(dades)
## Rows: 6,193
## Columns: 47
## $ I_D <dbl> 1, 2, 3, 5877, 8, 11, 9712, 12, 13,…
## $ FLAG_TSUNAMI <chr> "No", "Yes", "No", "Yes", "No", "No…
## $ YEAR <dbl> -2150, -2000, -2000, -1610, -1566, …
## $ MONTH <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DAY <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ HOUR <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MINUTE <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ SECOND <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ FOCAL_DEPTH <dbl> NA, NA, 18, NA, NA, NA, NA, NA, NA,…
## $ EQ_PRIMARY <dbl> 7.3, NA, 7.1, NA, NA, NA, NA, 6.5, …
## $ EQ_MAG_MW <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EQ_MAG_MS <dbl> NA, NA, 7.1, NA, NA, NA, NA, NA, NA…
## $ EQ_MAG_MB <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EQ_MAG_ML <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EQ_MAG_MFA <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EQ_MAG_UNK <dbl> 7.3, NA, NA, NA, NA, NA, NA, 6.5, 6…
## $ INTENSITY <dbl> NA, 10, 10, NA, 10, 10, NA, NA, NA,…
## $ COUNTRY <chr> "JORDAN", "SYRIA", "TURKMENISTAN", …
## $ STATE <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ LOCATION_NAME <chr> "JORDAN: BAB-A-DARAA,AL-KARAK", "S…
## $ LATITUDE <dbl> 31.100, 35.683, 38.000, 36.400, 31.…
## $ LONGITUDE <dbl> 35.50, 35.80, 58.20, 25.40, 35.30, …
## $ REGION_CODE <dbl> 140, 130, 40, 130, 140, 130, 140, 1…
## $ DEATHS <dbl> NA, NA, 1, NA, NA, NA, NA, NA, NA, …
## $ DEATHS_DESCRIPTION <dbl> NA, 3, 1, NA, NA, NA, NA, NA, NA, N…
## $ MISSING <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ MISSING_DESCRIPTION <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ INJURIES <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ INJURIES_DESCRIPTION <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DAMAGE_MILLIONS_DOLLARS <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ DAMAGE_DESCRIPTION <dbl> 3, NA, 1, NA, 3, NA, 3, 3, 3, 3, NA…
## $ HOUSES_DESTROYED <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ HOUSES_DESTROYED_DESCRIPTION <dbl> NA, NA, 1, NA, NA, NA, NA, NA, NA, …
## $ HOUSES_DAMAGED <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ HOUSES_DAMAGED_DESCRIPTION <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_DEATHS <dbl> NA, NA, 1, NA, NA, NA, NA, NA, NA, …
## $ TOTAL_DEATHS_DESCRIPTION <dbl> NA, 3, 1, 3, NA, NA, NA, NA, NA, NA…
## $ TOTAL_MISSING <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_MISSING_DESCRIPTION <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_INJURIES <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_INJURIES_DESCRIPTION <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_DAMAGE_MILLIONS_DOLLARS <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_DAMAGE_DESCRIPTION <dbl> NA, NA, 1, 3, NA, NA, 3, NA, NA, NA…
## $ TOTAL_HOUSES_DESTROYED <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_HOUSES_DESTROYED_DESCRIPTION <dbl> NA, NA, 1, NA, NA, NA, NA, NA, NA, …
## $ TOTAL_HOUSES_DAMAGED <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ TOTAL_HOUSES_DAMAGED_DESCRIPTION <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
Com podem veure, tenim 6193 registres i 47 variables. Aquestes dades, inclouen dades que es remonten a l’any 2150 BC fins a l’actualitat. Observem la quantitat de NA’s que hi ha!.
Fem una descripció d’aquestes variables, (https://medium.com/@dmukherjeetextiles/significant-earthquakes-2150-bc-to-2023-ad-e550c60fe863):
DIMENSIONS TEMPORALS
DIMENSIONS QUE DESCRIUEN LES CARACTERÍSTIQUES
DIMENSIONS GEOGRÀFIQUES
DIMENSIONS QUE INDIQUEN CONSEQÜÈNCIES HUMANES I MATERIALS
La resta són totals que de valors que ja s’han descrit anteriorment.
Les dades que farem servir un moment o altre són: YEAR, MONTH, DAY, HOUR, FLAG_TSUNAMI, FOCAL_LENGTH, EQ_PRIMARY, INTENSITY, COUNTRY, LOCATION_NAME, LATITUDE, LONGITUDE, DEATHS, DEATHS_DESCRIPTION, DAMAGE_MILLIONS_DOLLARS, DAMAGE_DESCRIPTION.
Podem fer una primera ullada i mirar com es distribueixen les dades recollides al llarg dels anys.
ggplot(dades,
aes(x = YEAR)
) +
geom_histogram(binwidth = 2,
color = "darkred",
fill = "darkred" )
Podem veure que en temps més recents tenim més dades. Això és raonable ja que a mida que tirem enrere tenim menys informació de les variables.
Si mirem els valors NA:
print('NA')
## [1] "NA"
colSums(is.na(dades))
## I_D FLAG_TSUNAMI
## 0 0
## YEAR MONTH
## 0 407
## DAY HOUR
## 561 2042
## MINUTE SECOND
## 2247 3364
## FOCAL_DEPTH EQ_PRIMARY
## 2965 1791
## EQ_MAG_MW EQ_MAG_MS
## 4872 3265
## EQ_MAG_MB EQ_MAG_ML
## 4391 6009
## EQ_MAG_MFA EQ_MAG_UNK
## 6179 5416
## INTENSITY COUNTRY
## 3378 0
## STATE LOCATION_NAME
## 5872 1
## LATITUDE LONGITUDE
## 54 50
## REGION_CODE DEATHS
## 1 4129
## DEATHS_DESCRIPTION MISSING
## 3649 6172
## MISSING_DESCRIPTION INJURIES
## 6172 4955
## INJURIES_DESCRIPTION DAMAGE_MILLIONS_DOLLARS
## 4769 5685
## DAMAGE_DESCRIPTION HOUSES_DESTROYED
## 1762 5411
## HOUSES_DESTROYED_DESCRIPTION HOUSES_DAMAGED
## 4495 5710
## HOUSES_DAMAGED_DESCRIPTION TOTAL_DEATHS
## 5265 4500
## TOTAL_DEATHS_DESCRIPTION TOTAL_MISSING
## 4164 6168
## TOTAL_MISSING_DESCRIPTION TOTAL_INJURIES
## 6167 4940
## TOTAL_INJURIES_DESCRIPTION TOTAL_DAMAGE_MILLIONS_DOLLARS
## 4760 5740
## TOTAL_DAMAGE_DESCRIPTION TOTAL_HOUSES_DESTROYED
## 2920 5380
## TOTAL_HOUSES_DESTROYED_DESCRIPTION TOTAL_HOUSES_DAMAGED
## 4416 5772
## TOTAL_HOUSES_DAMAGED_DESCRIPTION
## 5384
Podem veure que hi ha moltíssims valors NA, algunes variables d’entrada són inservibles, d’altres pràcticament el 50% dels seus valors són nuls, d’altres ronda el 20% i poques es podríen netejar directament. Només FLAG_TSUNAMI, YEAR i COUNTRY estan netes, (3 de 47 !). Amb aquest panorama, si féssim una neteja com cal ens quedaríem a les fosques. Per tant decidim netejar algunes variables i conviure amb NA’s d’algunes altres, ho veurem tot seguit.
Abans que res, eliminem les darreres 16 variables que no farem servir per res i estan plegades de valors nuls.
dades <- dades |>
select(1:31)
Després d’anar comprovant com disminueix el nombre de registres(files), fem el següent:
dades_1900 <- dades |>
filter(YEAR >= 1900)
dades_netes <- dades_1900 |>
filter(YEAR >= 1900 & !is.na(FOCAL_DEPTH) & !is.na(EQ_PRIMARY) & !is.na(HOUR) & !is.na(DEATHS) & !is.na(DAMAGE_DESCRIPTION))
Podem fer un resum de com ens ha quedat el nostre dataset, (ho farem
usant la llibreria skimr,
https://search.r-project.org/CRAN/refmans/skimr/html/skim.html).
library(skimr)
skim_without_charts(dades_netes)
| Name | dades_netes |
| Number of rows | 1135 |
| Number of columns | 31 |
| _______________________ | |
| Column type frequency: | |
| character | 4 |
| numeric | 27 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| FLAG_TSUNAMI | 0 | 1.00 | 2 | 3 | 0 | 2 | 0 |
| COUNTRY | 0 | 1.00 | 3 | 36 | 0 | 99 | 0 |
| STATE | 1108 | 0.02 | 2 | 2 | 0 | 8 | 0 |
| LOCATION_NAME | 0 | 1.00 | 4 | 52 | 0 | 994 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 |
|---|---|---|---|---|---|---|---|---|---|
| I_D | 0 | 1.00 | 6096.53 | 2291.86 | 2583.00 | 4569.00 | 5291.00 | 7674.50 | 10494.00 |
| YEAR | 0 | 1.00 | 1984.63 | 27.54 | 1900.00 | 1969.50 | 1990.00 | 2006.00 | 2020.00 |
| MONTH | 0 | 1.00 | 6.55 | 3.42 | 1.00 | 4.00 | 7.00 | 9.00 | 12.00 |
| DAY | 0 | 1.00 | 16.09 | 8.59 | 1.00 | 9.00 | 16.00 | 24.00 | 31.00 |
| HOUR | 0 | 1.00 | 11.32 | 7.11 | 0.00 | 5.00 | 11.00 | 18.00 | 23.00 |
| MINUTE | 3 | 1.00 | 28.92 | 17.42 | 0.00 | 14.00 | 28.00 | 43.00 | 59.00 |
| SECOND | 113 | 0.90 | 30.49 | 17.31 | 0.10 | 15.00 | 30.20 | 45.80 | 59.90 |
| FOCAL_DEPTH | 0 | 1.00 | 31.67 | 40.55 | 0.00 | 10.00 | 22.00 | 33.00 | 651.00 |
| EQ_PRIMARY | 0 | 1.00 | 6.38 | 0.97 | 2.10 | 5.70 | 6.40 | 7.10 | 9.50 |
| EQ_MAG_MW | 597 | 0.47 | 6.54 | 0.88 | 4.20 | 5.90 | 6.40 | 7.18 | 9.50 |
| EQ_MAG_MS | 356 | 0.69 | 6.48 | 0.97 | 3.50 | 5.80 | 6.50 | 7.30 | 8.80 |
| EQ_MAG_MB | 422 | 0.63 | 5.90 | 0.64 | 2.10 | 5.50 | 5.90 | 6.30 | 8.20 |
| EQ_MAG_ML | 1100 | 0.03 | 5.96 | 1.07 | 3.60 | 5.45 | 6.00 | 6.65 | 7.70 |
| EQ_MAG_MFA | 1134 | 0.00 | 6.30 | NA | 6.30 | 6.30 | 6.30 | 6.30 | 6.30 |
| EQ_MAG_UNK | 983 | 0.13 | 6.34 | 0.87 | 3.60 | 5.80 | 6.30 | 7.03 | 8.30 |
| INTENSITY | 636 | 0.44 | 7.89 | 1.85 | 2.00 | 7.00 | 8.00 | 9.00 | 12.00 |
| LATITUDE | 0 | 1.00 | 19.68 | 21.46 | -54.00 | 2.59 | 27.38 | 37.15 | 61.02 |
| LONGITUDE | 0 | 1.00 | 38.36 | 79.62 | -178.25 | -0.97 | 52.19 | 103.37 | 178.29 |
| REGION_CODE | 0 | 1.00 | 102.80 | 54.92 | 10.00 | 40.00 | 130.00 | 150.00 | 170.00 |
| DEATHS | 0 | 1.00 | 1875.49 | 15372.49 | 1.00 | 2.00 | 8.00 | 60.00 | 316000.00 |
| DEATHS_DESCRIPTION | 0 | 1.00 | 1.56 | 1.01 | 1.00 | 1.00 | 1.00 | 2.00 | 4.00 |
| MISSING | 1116 | 0.02 | 2410.89 | 9946.70 | 1.00 | 6.50 | 21.00 | 157.00 | 43476.00 |
| MISSING_DESCRIPTION | 1116 | 0.02 | 1.84 | 1.01 | 1.00 | 1.00 | 1.00 | 3.00 | 4.00 |
| INJURIES | 422 | 0.63 | 3676.10 | 34627.61 | 1.00 | 21.00 | 100.00 | 400.00 | 799000.00 |
| INJURIES_DESCRIPTION | 339 | 0.70 | 2.26 | 1.13 | 1.00 | 1.00 | 2.00 | 3.00 | 4.00 |
| DAMAGE_MILLIONS_DOLLARS | 758 | 0.33 | 1652.10 | 7791.53 | 0.10 | 5.00 | 35.00 | 400.00 | 100000.00 |
| DAMAGE_DESCRIPTION | 0 | 1.00 | 2.51 | 1.05 | 1.00 | 2.00 | 2.00 | 3.00 | 4.00 |
Les dades s’han reduït de 6193 a 1135 registres, (això és considerable i potser segons com pot donar lloc a dades esbiaxades, però pels nostres objectius considerem que no).
Nota: les dades que tenim del 1900 endavant
(dades_1900), sense netejar són 3720 registres, les farem
servir en algunes visualitzacions, així tindrem més punts
representatius, ara bé pel cas de l’estudi PCA necessitarem les dades
amb NA’s a zero ja que sinó tenim problemes.
Per les dades que es mostren no s’observen outliers, en el sentit de dades fora del seu rang, i tampoc tenim dades repetides. Podem visualitzar les distribucions de diferents variables:
library(cowplot)
anys <- ggplot(dades_netes,
aes(x = YEAR)) +
geom_bar(color = "black", fill = "cornflowerblue")
mesos <- ggplot(dades_netes, aes(x = MONTH)) +
geom_bar(color = "black", fill = "cornflowerblue")
plot_grid(anys, mesos)
dies <- ggplot(dades_netes,
aes(x = DAY)) +
geom_bar(color = "black", fill = "cornflowerblue")
hores <- ggplot(dades_netes, aes(x = HOUR)) +
geom_bar(color = "black", fill = "cornflowerblue")
plot_grid(dies, hores)
profunditat <- ggplot(dades_netes,
aes(x = FOCAL_DEPTH)) +
geom_bar(color = "black", fill = "cornflowerblue")
magnitud <- ggplot(dades_netes, aes(x = EQ_PRIMARY)) +
geom_bar(color = "black", fill = "cornflowerblue")
plot_grid(profunditat, magnitud)
morts <- ggplot(dades_netes,
aes(x = DEATHS)) +
geom_bar(color = "black", fill = "cornflowerblue")
danys <- ggplot(dades_netes, aes(x = DAMAGE_MILLIONS_DOLLARS)) +
geom_bar(color = "black", fill = "cornflowerblue")
plot_grid(morts, danys)
Podem veure com més cap a l’actualitat tenim més dades, no s’observa
cap mes, dia o hora preferents en que es produeixin més terratrèmols,
(es podria mirar de veure si el nombre de morts estan relacionats amb
l’hora del dia en que s’ha produït el sisme).
En quant a la profunditat i magnitud, podria semblar que tenim
outliers, però els valors estan a dins dels establerts i a part
recordem que hem fet una neteja exhaustiva.
La distribució de terratrèmols, tenint en compte la magnitud. (Ho
representem en un mapa, per fer-lo més gran fem servir:
https://ggplot2-book.org/maps.html, https://www.tidyverse.org/blog/2020/08/taking-control-of-plot-scaling/)
:
library(maps)
world <- map_data('world')
title <- "Mapa de Terratrèmols segons la magnitud"
wr <- map_data("world")
p <- ggplot(dades_1900) +
geom_map(data = world, map = world,
aes(map_id = region),
fill="white",
colour="#7f7f7f") +
expand_limits(x = world$long, y = world$lat) +
geom_point(aes(x = LONGITUDE, y = LATITUDE, color = EQ_PRIMARY),
size = 1) +
ggtitle(title) +
labs(x = "Longitud", y = "Latitud", color = "Magnitud") +
scale_colour_gradient(high = "#132B43", low = "lightblue")
pngfile <- fs::path(knitr::fig_path(), "mapa1.png")
library(ragg)
agg_png(pngfile, width = 25, height = 12, units = "cm", res = 300)
plot(p)
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
invisible(dev.off())
knitr::include_graphics(pngfile)
Aquest mapa mostra l’activitat sísmica més significativa. Es pot veure com els punts estan alineats indicant els límits de les plaques tectòniques.
Considerem ara la presència de tsunamis:
dades_1900$tsunami <- ifelse(dades_1900$FLAG_TSUNAMI %in% c("Yes"), "Yes", "No")
counts <- table(dades_1900$tsunami)
barplot(prop.table(counts),
col = c("cornflowerblue","darkred"),
main="Terratèmols i tsunamis",
legend.text=c("No Tsunami","Sí Tsunami"),
xlab ="Tsunamis", ylab = "Percentatge",
ylim=c(0,0.8))
S’observa que un 75% dels terratrèmols de les nostres dades no provoquen tsunamis, mentre que un 25% sí. Això pot indicar que una gran part dels sismes es produeixen a zones terrestres i l’altre restant en zones costeres o marítimes.
Això ho podem representar de nou en un mapa.
world <- map_data('world')
title <- "Mapa de Terratrèmols segons la presència de tsunamis"
wr <- map_data("world")
p <- ggplot(dades_1900) +
geom_map(data = world, map = world,
aes(map_id = region),
fill="white",
color="#7f7f7f") +
expand_limits(x = world$long, y = world$lat) +
geom_point(aes(x = LONGITUDE, y = LATITUDE, color = FLAG_TSUNAMI),
size = 1) +
ggtitle(title) +
labs(x = "Longitud", y = "Latitud", color = "Hi ha hagut tsunami?")
pngfile <- fs::path(knitr::fig_path(), "mapa2.png")
agg_png(pngfile, width = 25, height = 12, units = "cm", res = 300)
plot(p)
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
invisible(dev.off())
knitr::include_graphics(pngfile)
Podem mirar les morts provocades pels terratrèmols que produeixen tsunami o no.
death_ts <- dades_1900 |>
filter(FLAG_TSUNAMI == "Yes")
death_nts <- dades_1900 |>
filter(FLAG_TSUNAMI == "No")
par(mfcol = c(1, 2))
hist(death_ts$DEATHS, breaks = 50, main = "Morts amb tsunami", xlab = "Nombre de morts", col = "cornflowerblue")
hist(death_nts$DEATHS, breaks = 50, main = "Morts sense tsunami", xlab = "Nombre de morts", col = "cornflowerblue")
Com és d’esperar, segons els resultats anteriors, les morts per sismes sense tsunami són d’un ordre de magnitud més gran. S’observa també que tot i que els grans terratrèmols provoquen gran quantitat de morts, (tant si hi ha tsunami o no), sempre hi ha un mínim de morts acumulats a quantitats més petites.
Si mirem els morts considerant la magnitud del terratrèmol:
graph1 <- ggplot(dades_netes,
aes(x = EQ_PRIMARY, y = DEATHS)) +
geom_point(color = "darkred") +
labs(x = "Magnitud", y = "Nombre de morts")
graph2 <- ggplot(dades_netes,
aes(x = EQ_PRIMARY, y = DEATHS_DESCRIPTION)) +
geom_point(color = "darkred") +
labs(x = "Magnitud", y = "Grup de mortalitat")
# https://ggplot2-book.org/arranging-plots
library(patchwork)
graph1 + graph2
Podem veure, complementant el resultat anterior, que la mortalitat s’incrementa amb la magnitud del terratrèmol, tot i que aquest resultat també posa de manifest que la densitat de població és un factor que incideix en la mortalitat, això es pot veure amb la presència de sismes de gran magnitud i pocs morts. (Recordem que els grups de mortalitat són: 0 = Cap, 1 = Poques, 2 = Algunes, 3 = Moltes i 4 = Moltíssims).
Podem mirar la magnitud dels terratrèmols que provoquen tsunami i els que no en provoquen.
mag_tsu <- dades_netes |>
filter(FLAG_TSUNAMI == "Yes")
mag_no_tsu <- dades_netes |>
filter(FLAG_TSUNAMI == "No")
g1 <- ggplot(mag_tsu,
aes( x = EQ_PRIMARY)) +
geom_bar(color = "black", fill = "cornflowerblue") +
labs(title = "Magnitud amb tsunami",
x = "Magnitud")
g2 <- ggplot(mag_no_tsu,
aes( x = EQ_PRIMARY)) +
geom_bar(color = "black", fill = "cornflowerblue") +
labs(title = "Magnitud sense tsunami",
x = "Magnitud")
plot_grid(g1, g2)
Es veu que els terratrèmols que provoquen tsunami tenen magnituds cap a la banda alta, els que no provoquen tsunami hi tenim un ventall més ampli de magnituds, (tenint en compte amb les dades amb les que treballem).
Centrem-nos ara amb la relació entre magnitud i profunditat:
ggplot(dades_1900,
aes(x = FOCAL_DEPTH, y = EQ_PRIMARY, color = EQ_PRIMARY)) +
geom_point() +
labs(x = "Profunditat (km)", y = "Magnitud", color = "Magnitud")
S’observa la presència de terratrèmols de totes les magnituds per profunditats baixes (< 100 km), aquests són la majoria , però a partir dels 200 km de profunditat les magnituds són superiors a 6, i a partir dels 300 km la magnitud és superior a 7. (Caldria fer una mica de recerca per trobar alguna relació amb l’estructura de capes terrestres i límits de plaques tectòniques).
Considerem els terratrèmols més profunds (profunditat > 200 km) i situem-los al mapa.
earthq300 <- dades_1900 |>
filter(FOCAL_DEPTH >= 200)
world <- map_data('world')
title <- "Terratrèmols més profunds ( > 200 km)"
wr <- map_data("world")
p <- ggplot(earthq300) +
geom_map(data = world, map = world,
aes(map_id = region),
fill="white",
color="#7f7f7f") +
expand_limits(x = world$long, y = world$lat) +
geom_point(aes(x = LONGITUDE, y = LATITUDE, color = EQ_PRIMARY),
size = 3) +
ggtitle(title) +
labs(x = "Longitud", y = "Latitud", color = "Magnitud") +
scale_colour_gradient(high = "#132B43", low = "lightblue")
pngfile <- fs::path(knitr::fig_path(), "mapa4.png")
agg_png(pngfile, width = 25, height = 12, units = "cm", res = 300)
plot(p)
invisible(dev.off())
knitr::include_graphics(pngfile)
Hem anat fent visualitzacions en el mapa, però ara volem especificar la relació dels terratrèmols amb la geografia, és a dir països on n’hi ha més i a on són més greus.
Comencem amb amb els països on l’activitat sísmica és més gran, llistarem els 20 primers:
top20_quant <- dades_1900 |>
group_by(COUNTRY) |>
summarize(earthq_num = n()) |>
filter(COUNTRY != "") |>
arrange(desc(earthq_num)) |>
head(n = 20)
# https://sebastiansauer.github.io/ordering-bars/
ggplot(top20_quant,
aes(x = reorder(COUNTRY, -earthq_num), y = earthq_num,
fill = earthq_num)) +
geom_bar(position = "dodge", stat = "identity") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(title = "Països amb més nombre de terratrèmols",
x = "País", y = "Nombre de terratrèmols",
fill = "Nombre terratrèmols") +
scale_fill_gradient(high = "#132B43", low = "lightblue")
Volem ara llistar els 20 pitjors terratrèmols en magnitud, a quines zones han passat:
top20_mag <- dades_1900 |>
group_by(LOCATION_NAME , EQ_PRIMARY) |>
filter(LOCATION_NAME != "") |>
arrange(desc(EQ_PRIMARY)) |>
head(n = 20)
ggplot(top20_mag,
aes( x = reorder(LOCATION_NAME, -EQ_PRIMARY), y = EQ_PRIMARY,
fill = EQ_PRIMARY)) +
geom_bar(position = "dodge", stat = "identity") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(title = "Terratrèmols més forts",
x = "Zona", y = "Magnitud",
fill = "Magnitud") +
scale_fill_gradient(high = "#132B43", low = "lightblue")
Els terratrèmols que han causat més morts, (que no han de coincidir amb els anteriors, és a dir amb els de més magnitud).
top20_dead <- dades_1900 |>
group_by(LOCATION_NAME , DEATHS) |>
filter(LOCATION_NAME != "") |>
arrange(desc(DEATHS)) |>
head(n = 20)
ggplot(top20_dead,
aes( x = reorder(LOCATION_NAME, -DEATHS), y = DEATHS, fill = DEATHS)) +
geom_bar(position = "dodge", stat = "identity") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(title = "Terratrèmols més mortals",
x = "Zona", y = "Morts",
fill = "Morts") +
scale_fill_gradient(high = "#132B43", low = "lightblue")
I els més destructius, (danys materials):
top20_dam <- dades_1900 |>
group_by(LOCATION_NAME , DAMAGE_DESCRIPTION) |>
filter(LOCATION_NAME != "") |>
arrange(desc(DAMAGE_DESCRIPTION)) |>
head(n = 20)
ggplot(top20_dam,
aes( x = reorder(LOCATION_NAME, -DAMAGE_DESCRIPTION), y = DAMAGE_DESCRIPTION, fill = DAMAGE_DESCRIPTION)) +
geom_bar(position = "dodge", stat = "identity") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(title = "Terratrèmols més detructius",
x = "Zona", y = "Index de destrucció",
fill = "Index de destrucció") +
scale_fill_gradient(high = "#132B43", low = "lightblue")
Observacions:
D’aquestes darreres 4 gràfiques es poden observar algunes
característiques. Xina encapçala els països amb més nombre de
terratrèmols, Japó també. En quant a magnitud, destaquem Xile,
Indonèsia, Japó, Rússia (Kamtxatka) i Alaska. Si mirem el nombre de
morts, Haití és el primer de la llista i no apareix a les llistes
anteriors, això ens pot indicar que és una zona altament poblada i un
país molt pobre. En aquesta llista de més mortalitat, no hi apareixen ni
Kamtxatka ni Alaska (amb terratrèmols de gran magnitud), ja que són
zones amb una densitat de població molt baixa.
En quant al llistat de la destrucció, llevat del terratrèmol de San
Francisco a California a principis del segle passat, s’observa que en
general els països que hi apareixen són majoritàriament països amb
construccions més precàries, destacant però Japó que està més enrere,
que com és conegut, és un dels països amb les millors mesures
antisísmiques.
Podem completar aquestes visualitzacions la relació entre el grau de danys relacionat amb la magnitud.
# https://r-graph-gallery.com/2d-density-plot-with-ggplot2.html
# https://stackoverflow.com/questions/43515112/reversing-default-scale-gradient-ggplot2
top20_dam <- dades_1900 |>
group_by(DAMAGE_DESCRIPTION)
ggplot(dades_1900, aes(x=EQ_PRIMARY, y=DAMAGE_DESCRIPTION)) +
geom_bin2d(bins = 10) +
scale_fill_continuous(high = "royalblue4", low = "steelblue1") +
labs(x = "Magnitud" , y = "Grau de danys")
## Warning: Removed 1145 rows containing non-finite outside the scale range
## (`stat_bin2d()`).
D’aquesta visualització no podem treure gaires conclusions, ja hem vist que els danys i la magnitud depèn molt de la densitat de població i de la riquesa dels països, (primer món vs tercer món).
Creem una nova variable numèrica a partir de FLAG_TSUNAMI, que pendrà valor 1 o 0 segons sigui “Yes” o “No”
# Convertim Yes = 1, No = 0
dades_netes <- dades_netes |>
mutate(tsunami = ifelse(FLAG_TSUNAMI %in% c("No"), 0, 1))
Tot seguit cercarem correlacions en funció dels morts i altres
variables que poden estar relacionades, (inclosa tsunami
que acabem de crear).
# Utilitzem aquesta llibreria per fer servir la funcio multiplot()
library('Rmisc')
#Crearem primer una llista per mostrar les gràfiques de correlacions
#Crearem una llista per mostrar els atributs que interessen.
n = c("YEAR", "DAY", "HOUR", "tsunami", "FOCAL_DEPTH", "INTENSITY", "LATITUDE", "LONGITUDE", "EQ_PRIMARY", "DAMAGE_DESCRIPTION", "REGION_CODE")
dataAux <- dades_netes |>
select(all_of(n))
histList2<- vector('list', ncol(dataAux))
for(i in seq_along(dataAux)){
message(i)
histList2[[i]]<-local({
i<-i
col <- log(dataAux[[i]])
ggp <- ggplot(dataAux,
aes(x = dades_netes$DEATHS, y=col)) +
geom_point(color = "gray30") +
geom_smooth(method = lm, color = "firebrick") +
theme_bw() + xlab("Morts") + ylab(names(dataAux)[i])
})
}
multiplot(plotlist = histList2, cols = 4)
D’aquetes gràfiques podem veure:
La matriu de correlacions la podem visualitzar amb la funció
corrplot.
# https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html
library("corrplot")
n = c("YEAR", "DAY", "HOUR", "tsunami", "FOCAL_DEPTH", "INTENSITY", "EQ_PRIMARY", "LATITUDE", "LONGITUDE", "DEATHS", "REGION_CODE", "DAMAGE_DESCRIPTION")
factors <- dades_netes |>
select(all_of(n))
# https://stackoverflow.com/questions/14674431/error-in-eigencorr-infinite-or-missing-values-in-x-when-making-a-correlat
res <- cor(factors, use = "complete.obs")
corrplot(res,
method = "color",
tl.col = "black",
tl.srt = 30,
order = "AOE",
number.cex = 0.75,
sig.level = 0.01,
addCoef.col = "black")
No veiem cap correlació significativa, a part de les comentades anteriorment, que en pugui permetre eliminar cap de les variables.
No hem trobat cap correlació significativa de manera que poguem
eliminar variables, per tant considerem que el nostre model no té
informació que el pugui debilitar. Hem afegit una variable numèrica nova
tsunami que pren valors 0 /1, i és simplement la
codificació de la variable FLAG_TSUNAMI que pren valors
“Yes” / “No”. Les altres variables són exactament les mateixes que hem
mostrat després de la neteja.
Abans hem vist la relació entre magnitud i profunditat, el que farem ara és discretitzar la variable FOCAL_DEPTH, la profunditat. Ho farem seguint la divisió que fa USGS dels terratrèmols:
summary(dades_netes[,"FOCAL_DEPTH"])
## FOCAL_DEPTH
## Min. : 0.00
## 1st Qu.: 10.00
## Median : 22.00
## Mean : 31.67
## 3rd Qu.: 33.00
## Max. :651.00
dataAux["depth"] <- cut(dataAux$FOCAL_DEPTH, breaks = c(0, 70, 300, 700), labels = c("Shallow", "Intermediate", "Deep"))
head(dataAux$depth)
## [1] Shallow Shallow Shallow Shallow Shallow Shallow
## Levels: Shallow Intermediate Deep
plot(dataAux$depth,
main = "Nombre de terratrèmols per segment de profunditat",
xlab = "Segment de profunditat", ylab = "Quantitat",
col = "red4")
Tot seguit normalitzem (per la diferència), les variables per tal que cada variable contribueixi per igual al nostre anàlisi.
# Definim la funció de normalització
nor <- function(x) {(x -min(x))/(max(x)-min(x))}
# Guardem un nou dataset normalitzat
dades_netes$type<- NULL
n = c("YEAR", "DAY", "HOUR", "tsunami", "FOCAL_DEPTH", "EQ_PRIMARY", "LATITUDE", "LONGITUDE", "DEATHS", "REGION_CODE", "DAMAGE_DESCRIPTION")
dades_netes <- dades_netes |>
select(all_of(n))
dades_netes_nor <- as.data.frame(lapply(dades_netes, nor))
head(dades_netes_nor)
## YEAR DAY HOUR tsunami FOCAL_DEPTH EQ_PRIMARY LATITUDE
## 1 0.000000000 0.9333333 0.39130435 1 0.03840246 0.8513514 0.5651338
## 2 0.008333333 0.2666667 0.78260870 1 0.05069124 0.8243243 0.8224871
## 3 0.016666667 0.4000000 0.39130435 0 0.02304147 0.6486486 0.8233565
## 4 0.016666667 0.6000000 0.08695652 1 0.05069124 0.7297297 0.5912170
## 5 0.016666667 0.7000000 0.13043478 0 0.04608295 0.7567568 0.8162272
## 6 0.016666667 0.5000000 0.21739130 0 0.01382488 0.5810811 0.8242260
## LONGITUDE DEATHS REGION_CODE DAMAGE_DESCRIPTION
## 1 0.3148344 7.594961e-05 0.5000 0.6666667
## 2 0.8990557 5.379764e-05 0.1250 0.3333333
## 3 0.6362543 2.689882e-04 0.1875 1.0000000
## 4 0.2447166 6.325969e-03 0.5625 0.6666667
## 5 0.7136643 7.908253e-03 0.1875 1.0000000
## 6 0.7027259 1.543992e-02 0.1875 1.0000000
Nota: Aquest estudi el fem seguint la PAC1 i també STDA-PCA
Volem treballar amb noves característiques anomenades components principals, independents entre sí. Això ens permet representar el joc de dades en un nou sistema de coordenades, format per aquestes components principals. Aquest sistema està més ben adaptat a la distribució del joc de dades, de manera que recull millor la seva variabilitat, són noves variables que s’obtenen de combinacions lineals o barreges de les variables inicials, de manera que aquestes components principals, en particular les primeres retenen la major part de la informació.
Apliquem PCA a les columnes del nostre dataset, ho fem primer
executant la funció prcomp().
pca.earthq <- prcomp(dades_netes_nor)
summary(pca.earthq)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 0.4314 0.3715 0.3336 0.3093 0.2845 0.22201 0.19025
## Proportion of Variance 0.2511 0.1862 0.1502 0.1290 0.1092 0.06649 0.04883
## Cumulative Proportion 0.2511 0.4373 0.5874 0.7165 0.8256 0.89214 0.94097
## PC8 PC9 PC10 PC11
## Standard deviation 0.16748 0.10051 0.05801 0.04736
## Proportion of Variance 0.03784 0.01363 0.00454 0.00303
## Cumulative Proportion 0.97881 0.99243 0.99697 1.00000
La funció summary(), ens mostra la proporció de
variància aplicada al conjunt total de cada atribut. Podem veure que
l’atribut 1 explica el 0.2511 de variabilitat del total de dades, en
canvi l’atribut 7 explica només el 0.04883.
El següent histograma mostra el pes de cada atribut sobre el conjunt total de dades:
library('factoextra')
#Els valors propis corresponen a la quantitat de variació explicada per cada component principal (PC).
ev= get_eig(pca.earthq)
ev
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 0.186125657 25.1090145 25.10901
## Dim.2 0.138024978 18.6200615 43.72908
## Dim.3 0.111304141 15.0153253 58.74440
## Dim.4 0.095655537 12.9042728 71.64867
## Dim.5 0.080915617 10.9158051 82.56448
## Dim.6 0.049290469 6.6494599 89.21394
## Dim.7 0.036194481 4.8827645 94.09670
## Dim.8 0.028049490 3.7839761 97.88068
## Dim.9 0.010101961 1.3627905 99.24347
## Dim.10 0.003364748 0.4539165 99.69739
## Dim.11 0.002243182 0.3026132 100.00000
fviz_eig(pca.earthq)
Seguint el guió de la PAC1, volem utilitzar el mètode de Kàiser per decidir quina de les variables obtingudes farem servir. Aquest criteri manté totes les variables la variància de les quals és superior a 1.
# Calculem la variància dels components principals a partir de la desviació estàndard
var_earthq <- pca.earthq$sdev^2
var_earthq
## [1] 0.186125657 0.138024978 0.111304141 0.095655537 0.080915617 0.049290469
## [7] 0.036194481 0.028049490 0.010101961 0.003364748 0.002243182
Aquests valors obtinguts no ens deixen clar quines components hem de considerar. Per ajudar-nos a triar-les, escalem les dades, cosa que no havíem fet abans i tornem a calcular la variància.
# Escalem les dades
earthq_scale <- scale(dades_netes_nor)
# Calculem les components principals
pca.earthq_scale <- prcomp(earthq_scale)
# Mostramos la varianza de dichas variables:
var_earthq_scale <- pca.earthq_scale$sdev^2
head(var_earthq_scale)
## [1] 2.0645709 1.6218845 1.1125793 1.0508882 1.0385397 0.9398405
Després d’analitzar la variància i aplicant el criteri de Kàiser ens quedarem amb les components principals 1, 2, 3, 4 i 5 que són les superiors a 1. Aquest criteri té el problema de sobreestimar el nombre de factors, però malgrat això és el que aplicarem per analitzar els resultats.
Mostrem l’histograma de percentatge de variància explicat amb les dades escalades:
fviz_eig(pca.earthq_scale)
ev = get_eig(pca.earthq_scale)
ev
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.0645709 18.768827 18.76883
## Dim.2 1.6218845 14.744405 33.51323
## Dim.3 1.1125793 10.114357 43.62759
## Dim.4 1.0508882 9.553529 53.18112
## Dim.5 1.0385397 9.441270 62.62239
## Dim.6 0.9398405 8.544005 71.16639
## Dim.7 0.8861009 8.055463 79.22185
## Dim.8 0.7756779 7.051617 86.27347
## Dim.9 0.5961553 5.419594 91.69307
## Dim.10 0.5322558 4.838689 96.53175
## Dim.11 0.3815070 3.468245 100.00000
Aquesta taula s’interpreta de la següent manera, la suma de tots els valors propis és 11. Així considerant el primer valor propi, 2.0645709 i dividint-lo per 11, obtenim el 18.768827 percent de la tercera columna i la darrera és la suma dels valors anteriors acumulats.
Aquests valors propis es poden utilitzar per determinar el nombre de components principals a retenir després de la PCA (Kaiser 1961):
Observem que les 5 primeres components principals expliquen el 62.6% de la variació. És un valor una mica baix, però en aquest cas, per les dades que tenim ho donem per bo.
Així, després d’aplicar el mètode Kàiser ens quedem amb 5 components
principals. Per obtenir els resultats d’aquestes 5 variables del que ens
ha donat la sortida del PCA, fem servir la funció
get_pca_var() (factoextra package). Aquesta funció
ens dona una llista de matrius que contenen el resultat de les variables
considerades (coordenades, correlació entre variables i eixos, cosinus
al quadrat i contribucions).
var <- get_pca_var(pca.earthq_scale)
var
## Principal Component Analysis Results for variables
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the variables"
## 2 "$cor" "Correlations between variables and dimensions"
## 3 "$cos2" "Cos2 for the variables"
## 4 "$contrib" "contributions of the variables"
Els components de get_pca_var() es poden utilitzar en el diagrama de variables de la següent manera:
var.cos2 = var.coord \* var.coord.(var.cos2 \* 100) / (cos2 total del component).#Utilitzem els 5 components principals trobats abans
head(var$coord[,1:5],11)
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## YEAR -0.49726751 -0.26462329 0.22745564 0.25608947 -0.34108621
## DAY 0.07749664 0.02575068 -0.20154527 -0.67088882 -0.26413864
## HOUR 0.04217844 0.03367169 0.05360339 -0.07727443 -0.86173371
## tsunami 0.67205414 0.22971898 0.24013393 0.24746499 -0.01419688
## FOCAL_DEPTH 0.31095625 -0.16922316 0.40773756 -0.61222775 0.12059281
## EQ_PRIMARY 0.78900900 0.26567688 0.25134400 0.02580474 0.01789808
## LATITUDE -0.42314110 0.53895037 -0.35378868 -0.21624743 0.18312465
## LONGITUDE -0.36993791 0.52587765 0.53467986 0.03843890 0.01096036
## DEATHS 0.19837872 0.40833002 -0.18684736 0.14628881 -0.23133026
## REGION_CODE 0.44251807 -0.60327882 -0.35527994 0.13439894 -0.00372567
## DAMAGE_DESCRIPTION 0.29588280 0.54806669 -0.39758852 0.06936600 -0.08717138
La qualitat de representació de les variables en el mapa de factors s’anomena cos2 (cosinus quadrat, coordenades quadrades). Podem accedir al cos2 de la següent manera:
head(var$cos2[,1:5],11)
## Dim.1 Dim.2 Dim.3 Dim.4
## YEAR 0.247274973 0.0700254873 0.051736069 0.0655818182
## DAY 0.006005729 0.0006630974 0.040620496 0.4500918103
## HOUR 0.001779021 0.0011337825 0.002873324 0.0059713377
## tsunami 0.451656771 0.0527708094 0.057664303 0.0612389197
## FOCAL_DEPTH 0.096693789 0.0286364784 0.166249918 0.3748228133
## EQ_PRIMARY 0.622535208 0.0705842045 0.063173807 0.0006658845
## LATITUDE 0.179048389 0.2904675000 0.125166433 0.0467629512
## LONGITUDE 0.136854054 0.2765473012 0.285882554 0.0014775493
## DEATHS 0.039354115 0.1667334030 0.034911937 0.0214004159
## REGION_CODE 0.195822246 0.3639453363 0.126223834 0.0180630749
## DAMAGE_DESCRIPTION 0.087546630 0.3003771000 0.158076628 0.0048116422
## Dim.5
## YEAR 1.163398e-01
## DAY 6.976922e-02
## HOUR 7.425850e-01
## tsunami 2.015515e-04
## FOCAL_DEPTH 1.454263e-02
## EQ_PRIMARY 3.203412e-04
## LATITUDE 3.353464e-02
## LONGITUDE 1.201296e-04
## DEATHS 5.351369e-02
## REGION_CODE 1.388061e-05
## DAMAGE_DESCRIPTION 7.598849e-03
Podem visualitzar el cos2 de les variables en totes les
dimensions amb el paquet corrplot:
corrplot(var$cos2[,1:5], is.corr=FALSE)
També podem crear un diagrama de barres de variables
cos2 amb la funció fviz_cos2():
fviz_cos2(pca.earthq_scale, choice = "var", axes = 1:2)
Nota:
Un cos2 elevat indica una bona representació de la variable en el component principal. En aquest cas, la variable es dibuixada prop de la circumferència del cercle de correlació.
Un cos2 baix indica que la variable no està perfectament representada pels PC. En aquest cas, la variable està a prop del centre del cercle.
Per a una variable donada, la suma del cos2 de tots els components principals és igual a 1.
Si una variable està perfectament representada per només dos components principals (Dim.1 i Dim.2), la suma del cos2 en aquests dos PCs és igual a un. En aquest cas les variables es col·locaran en el cercle de correlacions.
Per a algunes de les variables, poden ser necessaris més de 2 components per representar perfectament les dades. En aquest cas les variables es situen dins del cercle de correlacions.
Resumint:
fviz_pca_var(pca.earthq_scale,
col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
Les contribucions de les variables en la comptabilització de la variabilitat d’un determinat component principal s’expressen en percentatge.
Les variables que estan correlacionades amb PC1 (és a dir, Dim.1) i PC2 (és a dir, Dim.2) són les més importants per explicar la variabilitat en el conjunt de dades.
Les variables que no estan correlacionades amb cap PC o amb les darreres dimensions són variables amb una contribució baixa i es poden eliminar per simplificar l’anàlisi global.
La contribució de les variables es pot extreure de la següent manera:
head(var$contrib[,1:5],11)
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## YEAR 11.97706360 4.31753848 4.6501017 6.24060839 11.202248904
## DAY 0.29089477 0.04088438 3.6510202 42.82965619 6.718011866
## HOUR 0.08616904 0.06990526 0.2582579 0.56821816 71.502800880
## tsunami 21.87654422 3.25367247 5.1829387 5.82734859 0.019407200
## FOCAL_DEPTH 4.68348110 1.76562994 14.9427477 35.66723912 1.400295553
## EQ_PRIMARY 30.15324882 4.35198712 5.6781397 0.06336397 0.030845348
## LATITUDE 8.67242616 17.90925926 11.2501134 4.44985018 3.229018205
## LONGITUDE 6.62869231 17.05098614 25.6954766 0.14060004 0.011567165
## DEATHS 1.90616436 10.28022667 3.1379279 2.03641220 5.152782374
## REGION_CODE 9.48488828 22.43965808 11.3451539 1.71883885 0.001336551
## DAMAGE_DESCRIPTION 4.24042734 18.52025221 14.2081223 0.45786432 0.731685955
Quant més gran sigui el valor de la contribució, més aportació de la variable hi haurà al component.
corrplot(var$contrib[,1:5], is.corr=FALSE)
Les variables més importants (que més contribueixen) es poden ressaltar a la gràfica de correlació de la següent manera:
fviz_pca_var(pca.earthq_scale, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
)
Les variables que apunten al mateix costat del gràfic estan
correlacionades de manera positiva, i les que apunten a costats oposats
del gràfic ho estan negativament.
Així, podem fer les següents observacions:
Hem treballat amb un dataset, Global Significant Earthquake Database from 2150BC, que és un llistat d’uns 6000 terratrèmols significatius que van del 2150 BC fins a l’actualitat. Aquests terratrèmols (registres) tenen un identificador únic i una sèrie de variables que el descriuen: variables temporals que en donen informació de quan van passar els sismes, variables que descriuen les característiques de terratrèmol, variables geogràfiques i variables que descriuen els efectes i danys sobre les persones i les infraestructures.
Al ser un dataset amb una gran presència de valors nuls, hem hagut de
descartar una gran quantitat de registres, de fet hem considerat la
sismicitat significativa a partir de l’any 1900, que és quan l’increment
de les dades disponibles va augmentant, i hem netejat les variables que
més ens han interessat, mantenint la resta. Això ha minvat la qualitat
de les dades, tot i així ens ha permès seguir amb l’estudi. (Val a dir
que amb paciència es podria cercar o construir un dataset més
complet).
Ens hem quedat al final amb 1135 registres, (tot i que per alguna
gràfica, en particular els mapes hem fet servir més punts per posar de
manifest algunes estructures o regularitats).
Hem vist que esl terratrèmols estan distribuïts seguint els límits de
les plaques tectòniques i que alguns d’ells, els que estan en zones
marítimes desencadenen tsunamis. Aquests provoquen un augment de morts,
però com que el nombre de terratrèmols sense tsunami és molt més
superior, la mortalitat també ho és.
S’ha vist que els terratrèmols amb tsunami generalment tenen magnitud
superior a 6 aproximadament.
També hem vist que l’increment de la magnitud del terratrèmol incrementa el nombre de víctimes, però a la vegada s’observa la presència de magnituds elevades amb pocs morts, això ens indica que la densitat de població pot ser un factor important.
Visualitzant la relació entre la magnitud i la profunditat hem pogut determinar que els terratrèmols a una profunditat superior a 300 km, la magnitud és més gran que 7.
Posant noms als països amb més nombre de terratrèmols i les zones amb els més forts, mortals i destructius hem deduït que factors com la densitat de població i riquesa i desenvolupament dels països afectats poden influir en el nombre de morts i danys ocasionats.
Finalment, amb l’estudi PCA hem vist una sèrie de relacions que hem comentat i també que amb una nova variable combinació d’altres amb una correlació inicial, es podria considerar com un índex de devastació o destrucció d’un terratrèmol (significatiu).